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In previous work, the authors documented the Multicomponent Ablation Thermochemistry 
(MAT) and Fully Implicit Ablation and Thermal response (FIAT) programs. In this work, 
key features from MAT and FIAT were combined to create the new Fully Implicit Ablation, 
Thermal response, and Chemistry (FIATC) program. FIATC is fully compatible with FIAT 
(version 2.5) but has expanded capabilities to compute the multispecies surface chemistry 
and ablation rate as part of the surface energy balance. This new methodology eliminates B ' 
tables , provides blown species fractions as a function of time, and enables calculations that 
would otherwise be impractical (e.g. 4+ dimensional tables) such as pyrolysis and ablation 
with kinetic rates or unequal diffusion coefficients. Equations and solution procedures are 
presented, then representative calculations of equilibrium and finite-rate ablation in flight 
and ground-test environments are discussed. 
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Nomenclature 


= reaction constant [SI units] 

= non-dimensional mass flux 
= specific heat, J/kg-K 
= atoms of element k in species n 
= vector of equation errors 
= reaction activation temperature, K 
= radiation view factor 
= recovery enthalpy, J/kg 
= static enthalpy, J/kg 
= identity matrix 
= number of gas-phase species 
= Jacobian matrix 

= 1 - K, number of non-base gas species 
= diffusion flux, kg/m 2 -s 
= equilibrium constant [pressure units] 

= number of elements 
= element or gaseous base species 
= number of condensed-phase species 
= molecular weight, kg/kmol 
= number of reactions 
= surface mass flux, kg/m 2 -s 
= molecular formula for a species 
= pressure, kg/rn-s 2 
= conductive heat flux, W/m 2 
= radiative heat flux, W/m 2 

= net reaction rate, kmol/m 2 -s 
= forward reaction rate [SI units] 

= temperature, K 
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environment temperature, K 
surface-normal velocity, m/s 
vector of unknowns 
surface mole fraction 
mass fraction 

diffusion-coefficient-weighted fraction 
absorptivity 

Kronecker delta (=1 if i = k, else =0) 
emissivity 

blowing reduction parameter 
density, kg/m 3 

heat transfer coefficient, kg/m 2 -s 
mass transfer coefficient, kg/m 2 -s 
Stefan-Boltzmann constant, W/m'-K 4 
reaction stoichiometric coefficient 
molecules of base species k in species n 
reaction temperature exponent 


subscripts 

c 

e 

f 

g 

i 

j 

k 

l 

m 


char 

boundary layer edge 
failing species 
pyrolysis gas 

any gaseous species (both j and k) 
non-base gaseous species 
element or gaseous base species 
condensed phase species 
reaction 
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n = all species (both i and l) superscripts 

w = surface value P = product species 

1 = unblown value R = reactant species 

T w = evaluated at wall temperature 

Introduction 

TN existing analysis codes for thermal protection system (TPS) materials, the ablation chemistry is modeled (at 
best) using tables that are created in advance for specified compositions of material and of atmospheric gas. In the 
past, thermochemical ablation tables for specific gas environments and thermal protection materials were generated 
using codes such as ACE or MAT. 1 ' 2 Each line of the table represents a minimal amount of information from a 
solution of the multispecies ablation chemistry with various applicable constraints. These tables are input and 
interpolated by TPS analysis codes such as CMA or FIAT 3 ' 4 that solve the time-dependent material response, but the 
detailed chemical information such as blown species concentrations is lost. 

This methodology is subject to errors and uncertainties from table interpolations. Additionally, for complex 
problems such as pyrolysis and ablation with kinetically-controlled rates or unequal diffusion coefficients, it is 
impractical to create and to interpolate tables with four or more independent variables. Furthermore, for the same 
materials, different tables are needed for each atmospheric or ground-test environment to be analyzed. This 
requirement is cumbersome and inconvenient. For example, the Orion and MSL Projects tested Phenolic 
Impregnated Carbon Ablator (PICA) in various arcjet environments, 5 and analyses were performed for vehicles and 
trajectories for both Earth and Mars. To perform analyses for this collection of test and flight environments, ablation 
tables were generated for 16 different mixtures of nitrogen, oxygen, argon, and carbon dioxide - which resulted in 
over 90,000 lines of input for the PICA material. 

In this work, we extensively revised both the MAT chemistry code and the FIAT ablation and thermal response 
code to enable the integration of the two codes. The multispecies chemistry is embedded within the surface energy 
balance (SEB) boundary condition of the new Fully Implicit Ablation, Thermal response, and Chemistry (FIATC) 
code. In FIATC the ablation and gas species blowing rates are calculated dynamically as part of the material 
response solution, without using B' tables.* Incorporation of this chemistry package within FIATC enables the 
calculation of more complex problems, as mentioned previously, and several future upgrades are planned including 
ablators with graded composition, pyrolysis gas chemistry and coking, and stacked pyrolyzing ablators. For the 
PICA example given above, the input is under 150 lines of data and independent of the gas composition. 

The purpose of this work is to document the equations and solution procedures in FIATC, then to present 
representative calculations of ablation in ground-test and flight environments. [In the following sections we define 
the general equations with material fail, unequal diffusion coefficients, and unequal heat and mass transfer 
coefficients. We have not yet implemented all of these models, but plan to do so for analysis of the Avcoat heatshield 
material for the Orion Crew Module .] 

I. Multicomponent Ablation Equations 

A general theory for calculation of the thermochemical ablation rate of pyrolyzing TPS materials was presented 
in Reference 2. This theory was incorporated into the MAT code that included capabilities for simultaneous 
ablation, pyrolysis, surface element constraints, non-equilibrium surface reactions, and material failure. Several 
changes to the equations and to the solution methodology were required for dynamic integration with the SEB 
boundary condition equation. This section describes in detail the governing equations and solution methodology for 
multispecies thermochemical ablation as implemented in FIATC. 

Figure 1 presents an illustration of the mass fluxes at the surface of a thermal protection material in a boundary- 
layer flow. A thin control volume is positioned just above the surface. The symbol m is mass flux, and the mass 
fractions of elements in the charred solid, the pyrolysis gas, and the boundary-layer edge gas are denoted as Ykc> 
V and Y ke , respectively. We allow one condensed species with mass fractions Y klf to fail. 

For steady conditions there is no accumulation of elements in the control volume; therefore, element 
conservation is 


■ Our implementation differs significantly from that of Reference 6 in which MAT was used as a subroutine to generate partial 
B 'tables at each time step within FIAT. 
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m g^kg + m c^kc J kw + (PV) W + m fYklf (1) 

which summed over all elements yields the mass balance 

m g + m c = (pV) w + m f (2) 

because the sum of the diffusion fluxes is zero. 

Following theory presented in References 1 and 3, the diffusion flux is approximated using a transfer potential 
methodology, specifically 


Jkw Pe u e^M (Z kw Zf, e ) (3) 

where Z is a diffusion-coefficient-weighted average of mass and mole fractions. For equal diffusion coefficients, Y 
and Z are equal. With this assumption, elimination of ( pV) w and j kw from Equation (1), division by p e u e C M , and 
rearrangement of terms yields the element flux balance 55 

(l + B’ g+ B' c - B’ f )Y kw = Y ke + B' g Y kg + B' c Y kc - B' f Y klf (4) 

where B' = m / p e ii e C M is a non-dimensional mass flux. The element flux balance insures that the correct elemental 
composition Y kw is obtained in the gas phase at the surface. 

We consider next the equations for the case of equilibrium chemistry in the control volume of Figure 1. We 
assume the system contains a mixture of K elements, / gaseous species, and L condensed (liquid or solid) species. 
Simple sums over the gaseous species provide the mixture pressure 

p-2 p i (5 ) 


the mixture enthalpy 


PMh w =2hiPi 

i 


( 6 ) 


and the gaseous mass fraction at the wall 


PMY kw = M k ^ j c ki P i 
i 


( 7 ) 


In these equations, the un-subscripted 2M is the mixture molecular weight. 

Next, we select K gaseous species (symbolically denoted as N k ) to represent base species for the elements in the 
system. Non-base gaseous species N j are represented by formation reactions of the base species 

^ v kjN k N j 

k 

where v k j is obtained by inverting the relation c k 'j = c k ' k v k j (see nomenclature).** These formation reactions have 
temperature-dependent equilibrium constants 

55 For unequal diffusion coefficients Equation (4) is more complex. Specifically, the right-hand side contains both Z and Y 
coefficients, and the left-hand side includes a sum over species pressures. Details, found in Reference 3, were omitted here to 
save a page of text. 

Base species must have nonzero partial pressures and possess a nonsingular stoichiometry submatrix c k ' k ■ There must be at 
least one gaseous species containing each element. Otherwise, the non-present element must be removed from the equation set. 
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X j^ = P jYl P k Vki ( 8 ) 

k 

that are calculated from the species thermodynamic data as a function of temperature T. 

Condensed species N / also are represented by fonnation reactions of the gaseous base species: 

^ j y u N k -N l 

k 

These formation reactions have temperature dependent equilibrium constants defined as 

Xm-xflV* ( 9 ) 

k 

The coefficients Xj account for the fact that L condensed species share the available surface area. If the mixture of 
condensed species is ideal, then Xj are surface mole fractions. With this interpretation, we add the constraint 

2 x i= l ( 1Q ) 

/ 

that allows condensed species to share the available surface area. 

In the preceding equations, the known quantities are: P, T, B' g , v kn , all element fractions Y, and thermodynamic 
data required for calculation of species enthalpies and equilibrium constants. The number of unknown variables and 
independent equations is summarized in Table 1. One additional equation is needed to determine the failure rate B'f . 

In the absence of material failure, we set B'f = 0 to close the equation set. If a failure temperature is specified for a 
condensed species, different scenarios are possible. If the fail species is the only surface species, then when T 
reaches 7f ai j a range of B'j is possible up to some maximum value that drives one to zero in Equation (4). Or if 

the failing species is one of multiple surface species, then a range of B'j-- is possible for any T — Tfail ■ In either case, 
it is possible to iterate using prescribed values of B'f until we find the maximum value that drives one Ykw to zero. 
The goal is to use the maximum failure rate when a range of failure rates is possible; otherwise, the solution is not 
unique. [This iteration is inefficient computationally. We are investigating what is the best way to implement 
material failure for carbon/glass systems.] 

MAT includes an option to specify M reversible or irreversible heterogeneous reactions in the form 

n n 


Table 1. Equations and unknowns 


Variables 

Number of 
such variables 

Equation Number 

Number of 
such equations 

Pi 

/ 

( 4 ) 

K 

x t 

L 

( 5 ) 

1 

Y kw 

K 

( 6 ) 

1 

M 

1 

( 7 ) 

K 

h w 

1 

( 8 ) 

I-K 

B'c 

1 

( 9 ) 

L 

B 'f 

1 

(10) 

1 

Total variables: 

I+K+L+4 

Total equations: 

I+K+L+3 
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Each reaction consumes the same condensed-phase species (typically carbon) denoted as Ny. The net forward 
reaction rate is expressed as 


R,„ = T- 


XmtX\r? L -X-'Ylpf 


The forward rate is written in Arrhenius form 


r m = A m T*"' exp (-E m IT) 


( 11 ) 


( 12 ) 


For irreversible reactions, the second tenn on the right-hand side of Equation (11) is omitted. There may be one or 
two gaseous reactants and products in each reaction. For carbon ablation a variety of reactions are possible. In this 
work we employ a modem model from Park: 6 

(-■(s') ■+■ o ^ CO 

C (S) +N^CN (13) 

3C (s) 


We omit the pseudo-reaction O2 + C — * 2 CO that is found in some older models. 

The chemical reactions are implemented as a conversion of one pseudo-element in the condensed species into 
another element of the same atomic weight in one of the product species. The quantity 


B 'kr = ' 


Mt 


Pe u e^ A 


■SSI 


ftmn 


p 

' M mn 


) c kn^-n 


which is the dimensionless consumption rate of element k by chemical reactions, appears as an extra term in the 
element flux balance. 


(l + B' g + B' c - B'f-)Y kw = Y ke + B' g Y kg + B' c Y kc - B'yY klf - B' kr (14) 

Examples of kinetic ablation of carbon phenolic will be presented in the Results section. 

II. Chemistry Solution Procedures 

The nonlinear equation set is solved by a Newton-Raphson procedure. There are advantages to using logs of 
some positive quantities as variables; therefore, the primary unknowns are chosen to be the vector 
X = [Y kw ,\nP k ,B' c ,lnPM,h w ,B'j-,lnPj ,lnX(]. The Appendix contains a summary of the equation errors and 
partial derivatives for the Newton-Raphson solution. The error vector E can be any convenient ordering of 
Equations (A1-A7). 

Given an initial guess for X, the error vector E and the Jacobian matrix J = t?E/o l X may be assembled. The 
Newton-Raphson procedure is iteration of the following: 

AX = -J -1 E (15a) 

X = X + AX (15b) 

In practice, the corrections to X are damped slightly, and provided J remains nonsingular, iteration of Equations 
(15ab) drives the error vector to zero, thereby yielding the solution X. Using the converged solution, partial 
derivatives with respect to any parameter may be determined as 
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c?X i <?E 

o' (parameter) r? (parameter) 


( 16 ) 


This feature is very useful, because at almost negligible additional computational cost, the derivatives of B' c and h w 
with respect to B' g and T may be calculated. 

For many ablative materials the total number of possible gas-phase species is large, especially in hydrocarbon 
mixtures from pyrolyzing ablators. As the number of gas-phase species increases, it becomes inefficient to factorize 
or to invert the large Jacobian matrix. However, most of the governing equations are the equilibrium-constant 
relations for non-base species. Each of these equations references only one non-base species (either Pj or A - ;); 
therefore, the Jacobian matrix contains an identity submatrix that results from the terms din OCj /o' In Py = 8 ^ and 
d\n%i IdlnPi' = 8 U ' in Equations (A6-A7). Therefore, these equations may be used to eliminate corrections InPj 
and \nX k that appear in the other equations. 

To implement this submatrix reduction, we define primary and secondary unknowns as listed in Table 2. The 
equations are ordered such that the equilibrium expressions are the final I + L-K equations. Equation (15a) is 
rewritten as 


'A 

B 

[AX ,1 


El' 

C 

I 

<N 

<l 

i 


E2 


(17) 


where A is a square matrix of dimension 2 K + 4 , and I is an identity matrix of size I + L-K . Submatrix reduction 
is the elimination of B to produce 


'A' 

O' 

AX/ 


K 

C 

I 

ax 2 


E2 


where A = A - BC and E'j = Ej - BE 2 . The solution is then 


AXj = -[A'] _1 Ej 
AX 2 =-E 2 -CAX! 


(IB) 


(19a) 

(19b) 


The advantage of the submatrix reduction is the large computational savings from inverting or factorizing the 
small matrix A 1 rather than the large matrix J, which more than offsets the expense of the simple algebra involved 
in Equations 17-19. Partial derivatives of solution variables (such as dB' c tdT ) may still be calculated, provided the 
vector of partial derivatives is carried through the algebraic operations described above. 

Convergence of the Newton-Raphson procedure will be adversely affected if, by poor choice of base species, 
some P k is very small compared with Pj for other species that contain element k. The solution of this potential 
problem is to dynamically reselect base species as required during the course of the solution. The approach used in 
FIATC is to select the species with the largest pressure containing each element. The matrix v kn depends on the 
choice of base species, but all other matrix operations are automated by use of pointer arrays for base and nonbase 
species. 

Trace species partial pressures or surface mole fractions may become too small to compute accurately. 
Provisions are made in FIATC to disregard species with mole fractions below a prescribed cutoff value (of 10'"). 
Non-base species at or below the cutoff value and with negative corrections to lnP ; or lnX; are removed from the 
equations and not considered in numerical convergence tests. 


Table 2. Variable sequencing for submatrix reduction 


Vector 

Designation 

Variables 

Length 

Xi 

primary 

Y kw ,in P k ,ln B' c ,ln PM,h w ,B'f 

2K + 4 

X 2 

secondary 

In Pj Jn Xj 

I + L-K 
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III. Surface Energy Balance 

We begin with a general form of the surface energy balance for unequal diffusion coefficients, that is slightly 
modified from that presented in Reference 3 (equation 78, page 48) 1 1 


p e u e C H (H r - h ew ) + p e u e C M | ^ (Z ie - Z iw )hJ™ - (B' c +B’ g - B' f )h w J 


+ m c h c + lilghg nt/hi F(76 w (f w l\c ) *7cond ^H’*?rad 


( 20 ) 


Here h ew is the enthalpy of the edge gas (frozen composition) evaluated at the wall temperature. The first term on 
the left-hand side is sensible convective heat flux; that is, the convective heat flux that would occur for a frozen 
boundary layer and a non-catalytic wall. The rest of the left-hand side includes all contributions from chemical 
reactions. The right-hand side terms are self-explanatory. 

In traditional CMA methodology, edge-gas tables provide h ew and ^ Z ie hJ w as functions of T and P, and 

thermochemical-ablation tables provide B h w , T, and ^ Z iw hJ w as functions of the independent variables P, B' g , 
and B' c . Both sets of tables apply to a specific edge gas composition. If the edge-gas composition changes as a 
function of time, the only feasible approach is to pick a fixed edge composition to use for the entire heating 
environment under consideration. Clearly this approach has questionable accuracy for trajectory -based 
environments. [ Varying edge-gas composition will be feasible with FIATC after we add the capability for unequal 
diffusion coefficients. The quantity h ew can be estimated using the chemistry subroutines provided the edge 
composition is specified as a function of time.] CMA allows input of a single value for C M /C H . For unequal 

diffusion coefficients. Reference 3 recommends C M = C H (Le) ~ where Le is the Lewis number. For equal 
diffusion coefficients ^Z le hJ w =h ew and ^ Z iw hJ w = h w ; however, the quantity h ew still is required. 

For equal diffusion coefficients and C M = C H , Equation (20) reduces to a more recognizable expression 

p e u e C H ^(H l . — h w ) + B c (h c — li w ) — Bf(h]f — h w )| + m g (h g — h w ) — Fos w (T w ~T CX) ') + q con d~ct w q ra d ( 21 ) 


For consistency with the coupling methodology described in the next section, the char and fail terms were moved 
inside the braces, whereas the pyrolysis gas terms were moved outside the braces. 

A blowing correction accounts for the reduction in heat transfer coefficient due to the injection of pyrolysis and 
ablation-product gases into the boundary layer. The blowing correction equation used by FIATC is 


C H ln(l + 2Afi') , n , (m c +m -m f ) 
= where B = — 

c hi 2A B’ p e ii e C M 


( 22 ) 


Here A is the blowing reduction parameter, and Ch /Cm is the ratio of the blown (ablating) to the unblown 
(nonablating) heat transfer coefficients. The fail mass flux does not contribute to the blowing correction. For laminar 
flow A = 0.5 or higher depending on the geometry and the ratio of molecular weights of the injected and boundary- 
layer-edge gas. For transitional or turbulent flow smaller values of A are used. 

IV. Coupling Implementation 


In FIATC the thermal-response solver calculates m g , T, g cond and density at all points within the solid 
including a surface node. The SEB from Equation (21) is a boundary condition. At each time step, the quantities 
H r , P, g rad , and p e u e C Hl are interpolated from environment tables. Then, the following steps are iterated. Based 
on current values, thermal properties for internal points and surface properties such as a w , e w , h c , h g , and h /f are 


In this formulation B'j is a portion of B' c ', that is, B' c represents the total solid ablation rate. However, in some codes B'j is 
treated as an independent quantity that is not included in B' c . It is important to determine which definition is used in any 
ablation tables with non-zero fail. 
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calculated. Using the current value of p e u e C H , the surface value of m g is converted into B' g . The values of B' g , P, 
and T w are input to the chemistry module that calculates B' c , B'p, and h w as well as derivatives with respect to B' g 
and T w . These quantities are loaded into the SEB equation of the thermal solver. The thermal solver uses a Newton- 
Raphson solution method, and the partial derivatives from the chemistry solution assist in convergence of the SEB 
equation. Typically 3 to 5 iterations are required to converge all equations including the SEB. 

For nonablating materials, the surface mass fluxes are zero, and the SEB simplifies to 

P ( ,U e C ij \ (H r — ll w ) = FO£ w (T w — Ty- ) + C] coru j — CL w Q r ad (23) 

In this case, the chemistry module calculates h w and dh w ldT w . 

Species and elemental information are required to implement the dynamic surface chemistry option in FIATC. 
This information is placed in the environment and material property input files. At the beginning of the environment 
data file, the user inputs one of the following nine keywords: Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune, 
Titan, or Elements. If one of the first eight keywords is input, FIATC uses the standard atmosphere for the selected 
planet or moon (Table 2). For Earth the table lists the composition of dry air,** and for Venus we omitted the trace 
amount of SCK If the keyword Elements is used, then FIATC reads the number of elements, their atomic weights, 
and their mass fractions in the edge gas. We provide a default database that contains JANNAF thermodynamic data 
of 91 gaseous species containing the elements H, He, N, O, Si, and/or Ar. This database is adequate for most flight 
or test atmospheres and for many thermal protection materials. 

In the material property file, each material has a type designation. Specifically, type 1 is reusable, type 3 is 
ablating but non-pyrolyzing, and type 4 is ablating and pyrolyzing. New choices in FIATC are types 5 and 6, which 
are the same as types 3 and 4, respectively, but without B' tables. For materials of types 5 and 6, new inputs are the 
number of elements, their atomic weights, and their mass fractions in the char and in the pyrolysis gas, followed by 
JANNAF data for the condensed-phase species under consideration, and finally data for up to five heterogeneous 
surface reactions. Calculations with dynamic surface chemistry may be performed with any materials of type 1, 5, or 
6. For surface materials of type 1, the simplified SEB from Equation (23) is used. [A material of type 2, pyrolyzing 
but non-ablating, also is feasible; however, we have not yet implemented this option with the dynamic surface 
chemistry. There has been little interest in this type of material.] 


Table 3. Known atmospheres 


Species 

Mole Fractions (from National Space Science Data Center fact sheets, at http://nssdc.gsfc.nasa.gov) 

Venus 

Earth 

Mars 

Jupiter 

Saturn 

Uranus 

Neptune 

Titan 

h 2 




0.898 

0.963 

0.825 

0.80 


He 




0.102 

0.0325 

0.152 

0.19 


co 2 

0.965 

0.00038 

0.9532 






ch 4 




0.003 

0.0045 

0.023 

0.015 

0.016 

n 2 

0.035 

0.78084 

0.027 





0.984 

0 2 


0.20946 

0.0013 






Ar 


0.00934 

0.016 






h 2 o 



0.00021 






CO 



0.0008 






nh 3 




0.00026 

0.000125 




NO 



0.0001 






Element 

Mass Fractions (normalized from the species tabulated above) 

Venus 

Earth 

Mars 

Jupiter 

Saturn 

Uranus 

Neptune 

Titan 

Hydrogen 



0.0000098 

0.8027572 

0.9133652 

0.6649572 

0.6401139 

0.0023185 

Helium 




0.1797725 

0.0606279 

0.2304166 

0.2909562 


Carbon 

0.2667604 

0.0001576 

0.2639229 

0.0158667 

0.0251909 

0.1046262 

0.0689299 

0.0069075 

Nitrogen 

0.0225654 

0.7551531 

0.0174533 

0.0016036 

0.0008160 



0.9907740 

Oxygen 

0.7106742 

0.2318084 

0.7038923 






Argon 


0.0128809 

0.0147217 







** For Earth the water vapor fraction is significant but highly variable in the troposphere. In the stratosphere (static pressure 
below about 10 kPa) the water fraction generally is below that of CCF. 
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V. Results 

[For the final paper we are considering a comparison of equilibrium and finite-rate PICA ablation for the 
following problems: 1) an arcjet environment, 2) thickness sizing for an Orion environment, 3) Stardust, and 4) an 
MSL environment. Avcoat calculations probably will be reported in a separate paper.] 

For this extended abstract, we consider two examples of ablation of phenolic impregnated carbon ablator 
(PICA). In previous work, an equilibrium ablation and thermal property model for PICA was developed. 8 The model 
was validated by comparison of predictions with arcjet data over a wide range of stagnation test conditions. 9 The 
predictions compared well with the data for surface temperature and recession for stagnation point heat flux above 
100 W/cm 2 . However, for heat fluxes below 80 W/cm 2 , corresponding to surface temperatures below 2000 K, the 
equilibrium model overpredicts the ablation rate, and a kinetic rate model for carbon ablation is recommended. 10 

In Reference 10, the lowest test condition was an air/argon mixture with measured stagnation heat flux and 
pressure of 40 W/cm2 and 4.4 kPa, respectively. The exposure duration was 200 seconds. The PICA model and 
facility calorimeter both had the iso-q shape, shown in Figure 2, with 10.16 cm diameter. This shape is used in 
testing in the NASA Ames arcjets, because the surface heating is relatively constant over the front face of the model. 
Following procedures described in Reference 9, facility data and flowfield analyses were combined to obtain 
boundary conditions for the material-response analysis. Specifically, in this case the centerline enthalpy and 
unblown heat transfer coefficient were estimated as 5.2 MJ/kg and 0.077 kg/m 2 -s, respectively. 

Reference 9 describes in detail some of the experimental uncertainties in arcjet testing. In order to assess the 
effect of these uncertainties on model predictions, our current procedure is to perform arcjet simulations using the 
nominal heating conditions and also using scaling factors of ±10% applied to the heat transfer coefficient. This 
procedure was followed using FIATC with the equilibrium ablation model for PICA, and the results are presented in 
Figure 3. The yellow-shaded region is the range of predicted recession. The black dots are the measured centerline 
recession for the two tested samples with error bars of ±0.5 mm. The equilibrium model overpredicts the measured 
recession by approximately 45%. 

The FIATC calculations were repeated using the reaction model for ablation of carbon phenolic in Equation (13) 
with rate parameters from Reference 7. The kinetic model predicts lower recession than the equilibrium model, and 
the range of predictions overlaps the data as shown in Figure 4. The predicted surface blowing rates are plotted in 
Figure 5. In this environment, the carbon ablation rate is predominantly a result of the reaction + O —■ * CO. The 

surface does not reach a steady state of ablation, although conditions are relatively constant over the last 50 seconds 
of heating. The char flux increases continuously, whereas the pyrolysis gas flux peaks within the first second of 
exposure. The total blowing rate peaks near 25 seconds, then decreases 5% over the next 175 seconds. 

The unblown heat transfer coefficient is considered to be constant for the entire test duration; however, the 
blown coefficient depends on the total surface mass flux according to the blowing correction in Equation (22). The 
maximum B ' g +B’ C is 0.174, which is results in a minimum blowing correction of 0.922; therefore, the minimum 
blown coefficient is 0.071 kg/m 2 -s (Figure 6). Also, except for the first 10 seconds of exposure, the blown 
coefficient is relatively constant. Therefore, for this specific arcjet case, approximate predictions could be obtained 
using kinetic ablation tables generated for a fixed blown coefficient of about 0.0713 kg/m 2 -s. However, in general, 
this approach is not accurate; instead, FIATC should be used to dynamically calculate the ablation rate without 
resorting to such assumptions. 

We compared FIATC predictions for kinetic and equilibrium ablation models with measured recession of PICA 
in arcjet tests over a broad range of stagnation conditions. The results are summarized in Figure 7. The equilibrium 
ablation model accurately predicts the recession for surface temperatures above 2000 K (or heat flux above 80 
W/cm 2 ). For surface temperatures below 2000 K, rate limitations to carbon oxidation become important, and the 
kinetic ablation model is required. Below about 1000 K the recession rate is negligible. 

Additionally, for pressure below about 7 kPa, the pyrolysis gas should be modeled as inert when it reaches the 
ablating surface. Pyrolysis gas is a multi-component hydrocarbon mixture that reacts with oxygen in the flowfield by 
a number of two- and three-body reactions. These homogeneous reactions have rates that are at least second in 
pressure. At sufficiently low pressure, these reactions become slower than the heterogeneous reaction for oxidation 
of solid carbon that is first-order in pressure, and it is incorrect to equilibrate the pyrolysis gas with the oxygen in the 
atmospheric gas. 

In our work to date, the kinetic model approaches, but does not reach, full equilibrium at high temperatures. 
Consequently, the kinetic model underpredicts the measured recession for conditions above 200 W/cm 2 . We are 
continuing our research in this area, because we would like to have a single model that may be used for the entire 
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range of flight conditions and which, of course, agrees with the entire set of stagnation test data. Our current 
procedure is to begin calculations using the kinetic model, then to switch to the equilibrium model for surface 
temperature above 2250 K. This hybrid model is inelegant, but consistent with the arcjet test data. 

The second example is a flight environment from a recent design cycle for the Orion Crew Module. Figure 8 
presents the margined heat fluxes at a critical location on the windward shoulder. The pressure history is shown in 
Figure 9. The vehicle first performs a skip through the upper atmosphere, then re-enters the atmosphere for a direct 
entry. Consequently there are two heat pulses. The first pulse, corresponding to the skip maneuver, has shorter 
duration but higher pressure and heat flux. The second pulse has longer duration but lower heat flux with negligible 
radiation. 

The total trajectory time is 1085 seconds; however, the ablative boundary conditions are turned off at 693 
seconds as the velocity drops below Mach 1 . After this time we apply a radiative cooling boundary condition that 
allows pyrolysis to continue during the cooldown portion of the trajectory. The TPS stackup is PICA bonded to a 
titanium-based structure using RTV-560 adhesive. As part of the analysis, FIATC calculates the minimum PICA 
thickness consistent with a maximum PICA-to-RTV temperature of 561 K. This temperature is the maximum reuse 
temperature of the adhesive. Orion flow environments are assumed to be fully turbulent; consequently, a blowing- 
reduction parameter A = 0.4 was used in the material response analyses 

FIATC calculations were performed using the equilibrium ablation model and the hybrid kinetic-equilibrium 
model. The predicted surface temperature (Figure 10) clearly shows the two heat pulses. For the hybrid ablation 
model, the equilibrium portion is required only during the hottest portion of the first heat pulse (the yellow-shaded 
region from 45 to 145 seconds), and the kinetic portion is used for the rest of the trajectory. The surface temperature 
is insensitive to the choice of ablation model, except at the tail end of the second heat pulse. After the second 
maximum in heat flux, the edge enthalpy decreases rapidly as the vehicle's velocity drops below Mach 1, and the 
atomic oxygen fraction also decreases. The kinetic model therefore predicts a low ablation rate, whereas the 
equilibrium model predicts steady recession at the diffusion-limited oxidation plateau ( B'c = 0.167). The wall 
enthalpy is lower for the latter case (owing to the presence of CO), and therefore the convective heating and surface 
temperature slightly larger. In both calculations the small kink at 693 seconds is a result of the switch in boundary 
conditions. 

The predicted ablation rate is shown in Figures 11. The kinetic model predicts lower recession rate in the cool 
period between the two heat pulses, and in the cooldown portion of the second heat pulse, as described previously. 
The 60% difference in predicted total recession (Figure 12) is remarkable. The equilibrium model predicts non-zero 
recession in the time interval between heat pulses; however, this inaccuracy is minor compared with its gross 
overprediction of the total recession in the second heat pulse. 

Interestingly, despite this large difference in predicted total recession, the calculated sizing to achieve the 
specified bond-line temperature was almost identical for the two cases: 5.71 cm for the equilbrium model and 5.76 
cm for the hybrid model. 

The surface blowing rates are presented in Figure 13. The char flux exceeds the pyrolysis gas flux during the 
highest heating portion of both heat pulses. The pyrolysis gas flux reaches a maximum in the first 50 seconds, then 
maintains a small value until beyond the end of the second heat pulse. There is no large increase in pyrolysis gas at 
the start of the second heat pulse, because any phenolic near the surface has already been fully charred. 

[Finally, Figure 14 shows the blowing rates of major species as a function of time...] 

[We will add a discussion of computational time for these calculations.] 

VI. Concluding Remarks 

[The final paper will include a conclusions section.] 
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Appendix: Error Equations and Derivatives for Ablation Chemistry 

The following is a summary of the error equations and partial derivatives used by the chemistry module within 
FIATC. 


Element flux balance 

E k = (l + B' g + B' c - B'f j Y kw - Y ke - B' g Y kg - B' c Y kc + B'jY kl/ - B' kr 
dE k / 3Y k ' w = 8 kk < (1 + B' g +B' c -B'f) 
dE k I dB' c = Y ^ - Y kc 
dE k !dB' f = Y klf -Y hv 

dE k IdlnPj = -dB' kr IdlnPj 
dE k Id In A - / = -dB' kr IdlnXj 
dE k ld\wT =-dB' kt . 1 3\nT 
dE k I dB' g = Jjtw - Y kg 

Wall mass fraction sum 


(Al) 


E k =M k ^c ki P i -PMY kw 

i 

dE k /3Y k ' w = -8 kk ’PM 
dE k Id In P t = M k c ki P; 
dE k !d\nPM = -PMY kw 


Pressure sum 


E = 2 P i~P 

i 

dE/dlnP^Pj 


(A2) 


(A3) 
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Mixture enthalpy sum 


E = ^ j h i P i -PMh w 

i 

dE / dlnPj = hjPj 
dE/dlnPM = -PMh w 
dE/dh w = -PM 
dE/dlnT = 2) CpiTPj 


Surface fraction sum 


E=^ x ,-\ 

I 

dE/d \nX l = X l 

Material failure rate [currently] 

E = Bf ~ (fif) prescribed 

dEj IdB'f = 1 

Equilibrium for non-base gaseous species 

E j = lnP J ~ E v kj ln Pk ~ In*, 

k 

dEj ld\nP k = -v k j 

dEj I d\n Pj < = d jj' 

dEj IdhvT = (2 v kj h k - hj)/RT 
k 


(A4) 


(A5) 


(A5) 


(A6) 


Equilibrium for condensed species 


E, = In A, v kl In P k - In*/ 

k 

dE t / d In P k =-v kl 

dE[ Id In X k = 6u’ 

dE l /dlnT=(2v k ,h k -h,)/RT 

k 


(A7) 
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Figure 2. Model shape for stagnation arcjet tests. TPS 

Figure 1. Mass fractions and compositions at the surface of samples have the same external shape as the calorimeter 
a thermal protection material in a boundary layer. used to measure the heat flux and stagnation pressure. 




Figure 3. Measured centerline recession and predictions of Figure 4. Measured centerline recession and predictions of 
equilibrium ablation model. kinetic ablation model. 
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Figure 5. Predicted surface blowing rates for nominal 
arcjet environment. 



Figure 6. Predicted blown heat transfer coefficient for 
arcjet environment. 
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Figure 7. Nominal ablation regimes for PICA. 


Figure 8. Margined heat fluxes at a critical shoulder 
location on the Orion Crew Module. 
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Figure 9. Pressure at a critical shoulder location on the 
Orion Crew Module. 




Figure 10. Surface temperature predicted by two ablation 
models for Orion environment. The equilibrium model is 
needed in the yellow-shaded region. 



Figure 11. Recession rate predicted by ablation models for Figure 12. Total recession predicted by ablation models for 
Orion environment. Orion environment. 
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Time, s 

Figure 13. Surface blowing rates predicted by hybrid 
ablation model. 


Plot to be inserted here 
species vs time 


Figure 14. Major blown species predicted by hybrid 
ablation model. 
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